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Proteins have been empirically linked to memory. If memory relates to protein structure, then each 
conformation would functionally code only one bit, making it difficult to explain large memories. 
Nor is there a simple way to relate memory to protein dynamics on current molecular dynamics 
(MD), which is memoryless. Here we point out that MD may be modified to involve memory 
ah initio without any new hypothesis: simply replace the electrostatic (Coulomb) force by the 
electrodynamic force — which is more accurate. We now need to solve functional differential equations 
(FDEs), instead of the ordinary differential equations (ODEs) currently solved in MD. Unlike ODEs, 
retarded FDEs are history-dependent: so memory is already present even at the level of interacting 
sites within molecules. The resulting increase in computational complexity is within the reach of 
current computers. While Amdahl's law does pose a challenge to parallelised time-stepping with 
this model, the compute-intensive part — the force calculation — may still be carried out in parallel. 
Thus, reformulating MD to use FDEs is feasible, and this could help to understand the possible 
dynamical basis of memory. 
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I. INTRODUCTION 

The functioning of biological macromolecules is linked 
to their structure and dynamics. The determination of 
structure through X-ray crystallography and NMR is ex- 
pensive, and these techniques provide no information 
about the dynamics of protein folding and protein inter- 
action. Molecular dynamics (MD) simulations are typi- 
cally used for this purpose and many computer pro- 
grams have been developed for MD, such as AMBER 

CHARMM GROMOS and NAMD 

Two key problems are (a) to determine the confor- 
mational structure (s) of protein molecules given the se- 
quence of amino acids, and (b) to determine the dy- 
namics of the folding process. The first problem is usu- 
ally studied by minimizing the free energy. Solving the 
Schrodinger equation for something as complex as the 
collection of particles constituting a protein molecule is 
too far beyond existing computational capabilities, so 
the second problem has long been studied by solving the 
Newtonian equations of motion for a prescribed force field 
involving van der Waals and electrostatic forces. 

Some features however seem hard to understand with 
this approach. First there is the classic Levinthal's para- 
dox: the conformational space to be explored is very vast 
in relation to the actual folding times that are observed. 
It has been hypothesised that this is because the energy 
landscape has a funneled character. While some attempts 
have been made to derive the funneled landscape from 



molecular dynamics, none of these is convincing. This 
hypothetical landscape is made more complex by possible 
multiple minima, evidence for which comes from prions. 

It is even harder to understand how long-term mem- 
ory formation in biological organisms links to protein 
molecules. Empirical studies have suggested that the for- 
mation of memory correlates with the formation of pro- 
tein molecules, and preventing protein molecules from 
forming impairs formation of long-term memories. But 
if a protein molecule can exist in only a few stable con- 
formations, then a large amount of functionally useful 
information cannot be stored in its structure, for each 
conformation corresponds to just one piece of functional 
information, irrespective of how that form was reached.^ 

Given the short folding times of protein molecules, the 
transitional state surely cannot be the seat of long-term 
memory. So, if protein molecules store memory struc- 
turally, then multibit information would require multi- 
ple conformations — to store memory on the scale of bio- 
logical organisms, a very large number of conformations 
would need to be associated with each protein molecule. 
If protein molecules do indeed admit such a large num- 
ber of stable conformations, then we are back to a situa- 
tion reminiscent of Levinthal's paradox, and we need to 
suppose that the energy landscape has a funnel with nu- 
merous minima or perhaps a multi-funnel landscape f^. 
Instead of thus accumulating hypotheses it seems better 
to look for a simpler solution. 

On the other hand, it is not immediately clear how 
memory can be present in the dynamics of protein in- 
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teractions either. Thus, MD simulations currently solve 
Newtonian equations of motion which are ordinary dif- 
ferential equations (ODEs). The solution of a system 
of ODEs is uniquely determined by the initial data (or 
data at any one instant of time), irrespective of past 
history. Thus, the time evolution of a system modeled 
by ODEs must be memoryless or history-independent. 
Hence, current MD, which involves memoryless dynam- 
ics, cannot readily explain memory-dependent molecular 
interactions in a simple way. It is possible, of course, 
that, as in artificial neural networks, memory is holistic 
and intrinsically beyond the reach of MD which is reduc- 
tionist. 

Before concluding that a simple explanation for mem- 
ory is intrinsically beyond the reach of molecular dynam- 
ics, it seems worthwhile to explore a more sophisticated 
mathematical model for MD which admits memory ab 
initio. The suggestion to explore such a model is not in- 
tended to exclude more holistic models of memory, but 
to complement them. Now, given the vast amount of 
experimental data that is currently available, it is not 
difficult to invent a variety of ad hoc mathematical mod- 
els which will fit some of it. However, that would go 
against the spirit of MD which has the great theoretical 
virtue that it proceeds directly from established physics, 
without invoking fresh hypotheses, except as simplifying 
assumptions. 

The modified model for MD, that we propose, needs no 
additional hypothesis. If existing physics is correctly ap- 
plied, the electrodynamic interaction between two mov- 
ing electrical charges is history dependent, hence already 
involves memory. If this feature is incorporated into 
molecular dynamics at the outset, it may help to un- 
derstand how protein molecules relate to memory. 

The argument, in outline, is very simple. Currently, 
the electrostatic (or Coulomb) potential is used in MD 
for long-range forces. However, the interacting sites are 
usually in motion so the Coulomb force is only an ap- 
proximation to the full electrodynamic force. If one uses 
the full electrodynamic force, the equations to be solved 
are functional differential equations (FDEs) Unlike 
ODEs the solution of these FDEs cannot be uniquely de- 
termined merely by initial data. Assuming these FDEs 
to be retarded, one must prescribe some part of the past 
history [1, [13, EH- So some rudimentary memory is al- 
ready present at the level of the constituents. 

In the case of two interacting elementary charges, a 
proton and electron, in a hydrogen atom, say, the system 
memory or the exact time interval over which past data is 
required is tiny, being of the order of a deci femto second. 
However, even in this case the difference between FDEs 
and ODEs is significant since solutions of retarded FDEs 
may exhibit qualitative behaviour that is impossible for 
solutions of ODEs. 

Further, the system memory, or the time period over 
which the past history of the system must be specified, 
increases with the number of interacting particles and the 
scale of the system. Thus, we have a simple reason to ex- 



pect more complex molecules and collections of molecules 
to exhibit longer-term memory. With currently available 
computational resources, computing solutions of FDEs 
for such large collections may be contemplated. 
This argument is elaborated below. 



II. EXISTING AND PROPOSED MOLECULAR 
DYNAMICS 



A wide variety of force fields have been used in molecular 
dynamic simulations of biological macromolecules [l2| . 
These force fields are typically derived from a potential 
represented as a sum over bonded and non-bonded pairs 
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The interaction potential V^j between non-bonded 
pairs is represented as a sum of a short-range van der 
Waals potential (usually a 12-6 Lennard-Jones potential) 
and the long-range electrostatic (Coulomb) potential: 



(1) 



Here, A and C are constants, qi and qj are charges, while 
Tij is the distance between them. 

Our concern here is with the second (electrostatic) 
term, on the right hand side of ([l]). No one earlier seems 
to have questioned the validity of using the Coulomb 
potential — which is used across various force fields in a 
variety of cornputer programs like amber CHARMM 
GROMOS Q, NAMD Q, etc. Currently used for MD 
simulations. 

However, it is elementary that the Coulomb poten- 
tial applies only to static charges, whereas in molecu- 
lar dynamics the interacting charges are in relative mo- 
tion. In the electrodynamic case, with moving charges, 
the force Fy on charge i due to charge j is actually ob- 
tained not from the Coulomb potential but from the (re- 
tarded) Lienard-Wiechert potential and the Heaviside- 
Lorentz force law (usually called just the Lorentz law), 
to give [11: 



R 



47reo (R • u)3 



{[(c^ - v^)u + Rx {u X a)] (2) 



R X [(c^ - v'^)u + R X (u X a)] | 



Here, charge qi is located at Vi{t) at time t, while the 
position of the other charge qj at time t is given by Tj (t) , 
R = ri(t) — rj(tr), and R is its norm. In the preceding 
expression, tr is the retarded time (the time at which the 
backward null cone with vertex at ri{t) meets the world 
line Tj {t) of the other charge) , and is given implicitly by 
the equation ||ri(t)— rj(tr)|| = c{t—tr), with ||'|| denoting 
the 3- vector norm, and c being the speed of light. Further 
(with dots denoting time derivatives), — ri{t), u — 
V = cR — V, and it is understood that v = rj{tr) 
Yj{tr) are the velocity and acceleration of the 
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charge qj at the retarded time tr. A similar expression 
gives the force Fji exerted on the charge qj by the charge 
qi. UnUke the electrostatic case, F^-j ^ F^ , in general. 
(Relativistic velocity effects are here ignored as irrelevant 
to molecular dynamics, although the conclusions apply a 
fortiori to the relativistic case.) 

The approximation of the full electrodynamic force ^ 
based on the Lienard-Wiechert potentials, by the sim- 
pler electrostatic force based on the Coulomb potential 
([I]) has a long history, going back to the days of the 
Rutherford model, and the Bohr atom. The general 
expectation among physicists has been that the use of 
the electrostatic force greatly simplifies calculations, and 
would not upset any key qualitative conclusion. Though 
long-standing, and widespread, this expectation is math- 
ematically incorrect, and has never been actually tested. 
Testing would require that one obtain the actual solu- 
tion of at least the 2-particle equations with the full 
electrodynamic force and check whether it is ap- 
proximately the same as the solution obtained with the 
electrostatic inverse square law force. Prior to the ad- 
vent of high-speed computers, doing this was well-nigh 
impossible: the one-particle case [15| being complicated 
enough, despite various attempts in the previous century 
0, Ea, H M, H mi, m, no actual solutions of the 2- 
particle equations with the full electrodynamic force were 
published, except [23| in oversimplified situations of little 
practical interest. 

The use of the full electrodynamic force ^ makes a 
fundamental mathematical difference, for it leads to the 
formulation of the n-body problem {n > 2) as a system 
of FDEs. 



III. FDES VS OTHER METHODS 

To see how FDEs differ fundamentally from ODEs, 
consider, for example, the simple FDE 
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It is easy to verify that both cos t and sin t are solutions 
of ([3]). Since the equation is linear, a cost + bsint is 
a solution of ^ for arbitrary constants a, b, and it is 
clear that the values of both constants a and b cannot 
be determined by a single initial condition a;(0) = a:o; 
say. The situation is shown in Fig. [TJ there is an infinity 
of non-unique solutions if we prescribe the state of the 
system at only one instant of time. To obtain a unique 
solution one must specify the past history of the system. 
A system modeled by FDEs, in effect, has memory. 

For the interaction of a proton and an electron within 
the classical hydrogen atom, the time interval over which 
past history must be prescribed is quite small. However, 
FDE exhibit qualitatively novel features, such as time 
asymmetry, and breakdown of phase flow as illustrated 
in Fig. [21 These novel features are, in principle, math- 
ematically impossible for solutions of ODEs. Therefore, 
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FIG. 1: Insufficiency of initial data. The figure shows 
three different solutions of a retarded FDE x{t) — x{t — 1), 
corresponding to three different past histories (prescribed for 
t < 0). All solutions have the same initial data {x{0)). The 
solution of a retarded FDE cannot be determined uniquely 
merely by prescribing initial data. 



regardless of the smallness of the time interval over which 
past history must be prescribed, we can expect qualita- 
tively new features to emerge from the application of this 
new model to molecular dynamics. 

It is possible to conceptualise the electrodynamic n 
body problem in an alternative way, using flelds. Since 
this has often been done in the past, we clarify how 
the field picture relates to our approach. To determine 
the field acting on a given particle, we must determine 
the total field generated by all other particles. To do 
this, it is necessary to solve Maxwell's equations. To 
solve these partial differential equations (FDEs), it is 
necessary to provide initial (Cauchy) data by prescrib- 
ing the electromagnetic field over all space at one instant 
of time (i.e., over a Cauchy hypersurface) . Assuming re- 
tarded Lienard-Wiechert potentials, the electromagnetic 
field due to a system of n-particles, at any instant of time, 
depends upon the past motions of those particles [23 |. 
Hence, prescribing this Cauchy data requires a knowl- 
edge of the entire past history of the positions, velocities. 
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FIG. 2: History Dependence of FDEs. The figure shows 
three different solutions of a retarded FDE x{t) — b{t)x{t — Vj, 
for a suitable choice of the function 6(t), as described in 
Three different past histories prescribed for t < G lead to 
three different solutions all of which coincide for f > 1. Such 
a phase collapse is impossible with ODEs where trajectories 
in phase space can never intersect. Because of this phase 
collapse, retarded FDEs, unlike ODEs, cannot, in general, be 
solved backward, from prescribed future data. 



and accelerations of the n particles producing the fields 
(Fig. 3). 

The field picture, therefore, only hides the dependence 
on the past, made expHcit in the particle picture. In the 
field approach to the electrodynamic n-body problem, we 
are required to solve a coupled system of ordinary and 
partial differential equations (PDEs). To determine the 
motion of one particle, we need to solve the Newtonian 
ODEs with the Heaviside-Lorentz force due to the fields 
generated by other particles. Those fields are determined 
from particle motions using Maxwell's equations. Com- 
pared to this coupled system of ODEs + PDEs in the 
field picture, it is, currently, computationally more con- 
venient to solve the FDEs of the particle picture for rea- 
sons already discussed elsewhere in detail in Which 
computational technique is used to solve the equations is 
not, of course, relevant to the model of memory that is 
being proposed here. 



IV. THE DIFFERENCE 

To ascertain the exact difference made by the use of 
the full electrodynamic force in place of the electrostatic 
force, we have computed [131 the first numerical solutions 
of the retarded FDEs of the 2-body problem, in a realis- 
tic physical context, using the full electrodynamic force, 
but neglecting radiation damping. In the case of the 
classical hydrogen atom, the solution with the full elec- 
trodynamic force (but without radiative damping) differs 
from the Coulomb orbit, prescribed as past data. This 
is summarily shown in Fig. IH The difference relates to 
an unexpected 'delay torque' [131 which arises because 
the full electrodynamic force depends upon the past mo- 
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FIG. 3: Relation of FDE method to ODE-|-PDE 
method. Assuming only two particles and retarded propaga- 
tors, the electromagnetic fields at any point a; on a Cauchy hy- 
persurface relate to past particle motions at points Xa and xt, 
where the backward null cone with vertex at x meets the world 
lines of the two particles a and b. As x runs over the hyper- 
surface, the points Xa and Xb will, in general, cover the entire 
past world-lines of the two particles. Thus, for solutions of the 
2-body problem, the field picture and the PDE+ODE method 
requests more information about the past than is practically 
needed by the particle picture and the FDE method. 

tion of the other particle. It is impossible for any central 
force (like the electrostatic force) to have such a torque. 
Thus, the electrodynamic interaction of two charged par- 
ticles involves complexities that cannot be captured by 
the simple Coulomb potential. 
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FIG. 4: Difference between FDE and ODE solution for 
the 2-particle problem. For the case of the classical hydro- 
gen atom, the figure plots the time evolution of the difference 
between the electron orbit with the full electrodynamic force 
and its orbit on the Coulomb force. The zero difference part 
corresponds to the past history, prescribed using the Coulomb 
orbit. 

This suggests that a far richer variety of behaviour 
can be modeled if the MD force field is reformulated by 
replacing the electrostatic force in (1) by the full electro- 
dynamic force (2). 

What are the computational costs of solving FDEs, in- 
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stead of the usual ODEs? Stiffness considerations can be 
set aside since they apply equally to FDEs and ODEs. 
The relative increase in complexity comes from the cal- 
culation of the electrodynamic force (2) instead of the 
Coulomb force in (1). However, once the retarded time 
tr is known, this only adds a constant number of floating 
point operations per pair of interacting sites. 

Thus, the basic increase in complexity comes from hav- 
ing to solve a nonlinear equation to determine tr for each 
(ordered) pair. This is not as bad as it sounds, since, 
between time steps, each value of tr would be expected 
to change by only a small amount. So, with the previ- 
ous value of tr as the starting guess, we can expect quick 
convergence in one or two steps. Thus, for practical pur- 
poses, the force calculation would asymptoticahy remain 
at most O(n^), though with a different constant. Further 
reduction of complexity by means of a cutoff is considered 
below. 

Memory requirements would also increase, since some 
part of the past history of each interacting particle / site 
would have to be kept in memory to avoid excessive swap- 
ping and interpolation. Exactly how much of the past 
history needs to be retained in memory depends upon the 
specific algorithm used to solve the FDEs, and whether 
or not it permits step sizes larger than the interval of 
retardation. 

On the whole, it is reasonable to expect that, for most 
existing MD computer programs, the resulting increase 
in computational complexity can be handled with exist- 
ing computers — an exceptional case is namd2 [5i], and 
the parallel versions of AMBER. While the solution of 
ODEs may be parallelised with reasonable efficiency [2H |. 
the history dependence of FDEs may be expected, by 
Amdahl's law, to pose a serious chahenge to parallel 
computing. However, Amdahl's law restricts only par- 
allelised time-stepping. The force calculation, which is 



the compute-intensive part, can still be done in parallel. 

The complexity can be reduced, as in the Coulomb 
case, by applying a long-range cutoff. This would reduce 
the memory requirements as well, though it is not so clear 
in the present context that this is necessarily desirable. 
The range of the fuh force, however, is larger, especiahy 
if we take into account the radiation damping. Until 
now, it was impossible to take into account the effects of 
radiation damping in a many-body problem, due to vari- 
ous long-standing difficulties like preacceleration [l5l.[lB|. 
However, these difficulties were recently addressed in a 
satisfactory way, using FDEs [2^1 • Theoretically, or with 
a long-term outlook, it is a significant advantage of this 
proposal that the effects of radiation damping can also 
be included, if so desired. Inclusion of radiation damping 
would, however, add to the complexity by introducing a 
new source of stiffness in the problem. 

Finally, we note that FDEs have long been linked to 
quantum mechanics, on the structured-time interpreta- 
tion of quantum mechanics [131 • A more detailed exami- 
nation of exactly how the FDE approach relates to quan- 
tum mechanics in the context of MD will be considered 
in a subsequent paper. 



V. CONCLUSIONS 

It is desirable to reformulate MD to use the fuh electro- 
dynamic force, which is more accurate than the Coulomb 
force. This involves solving FDEs, which is currently 
computationally feasible. This reformulation allows an 
exploration of a qualitatively richer set of ways in which 
biological macromolecules can interact at long range, and 
may help to understand the possible dynamical origin of 
biological memory. 
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